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Abstract 

The rough Heston model has recently attracted the attention of many finance practitioners and researchers as 
it maintains the basic structure of the classical Heston model while having descriptive capabilities in terms 
of micro-structural foundations of the market. Using the fact that the characteristic function of log-price in 
this model could be expressed in terms of the solution of a nonlinear parametric fractional Riccati differ¬ 
ential equation not admitting a closed-form solution, devising efficient numerical schemes for pricing and 
calibration under this model has become a crucial need in the computational finance community. Although 
the fractional Adams method has been used in most of the recent studies on the rough Heston model, this 
method suffers from some stability and convergence issues in treating the problem. In this paper, we present 
a numerical method based on Newton-Kantorovich quasi-linearization to solve the nonlinearity issue fol¬ 
lowed by spectral collocation based on “poly-fractonomials” to approximate the fractional derivative in an 
accurate and efficient manner. We provide sufficient conditions under which our method is convergent and 
the order of convergence is also obtained. In order to guarantee the specified convergence rate, we first prove 
some regularity results on the linearized problem and then employ the proposed scheme to solve a practical 
calibration problem from the SPX options market. The efficiency of the proposed method is illustrated by 
comparing the obtained results with those of the fractional Adams method as well as a fast hybrid scheme 
based on fractional power series expansion. 

Keywords: Rough Heston model. Fractional nonlinear Riccati differential equation. Spectral collocation, 
Newton-Kantorovich quasi-linearization. 

2010 MSC: 91B25, 91G60, 91G20, 34A08, 34A34 


1. Introduction 

The recent decade has seen a phenomenal growth in the number and variety of stochastic models for 
price fluctuations of financial assets. Most of these models are able to capture features like time-varying 
volatility, dynamic fat-tailed distributions and leverage effect very efficiently. Among them, those with 
explicit characteristic functions for the log-price relative are of special interest to finance practitioners and 
academia as they will result in efficient option pricing and calibration procedures [9, 46]. 

Continuous-time stochastic volatility (SV) models are among the flexible and widely used parametric 
asset price models in which the volatility process follows a stochastic differential equation driven by Brow¬ 
nian motion (see e.g. [25, 28, 43] for a review). However, empirical studies of market data show that both 
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historical and implied volatility time series display characteristics reminiscent of a fractional Brownian mo¬ 
tion (fBm) rather than a random walk process (see e.g. [11]). In this regal'd, the family of rough fractional 
stochastic volatility (RFSV) models [29, 20] were introduced into the field with the aim of overcoming the 
limitations of the classical SV models in describing the term structure of volatility smiles [10, 11], 

The use of fractional calculus in financial and economic modeling is not that new (see e.g. [3] for a 
comprehensive overview on the application of fractional tools in modeling long-memory processes, asset 
return distributions and option prices in the econometrics literature). The fractional stochastic volatility 
models were proposed in the late 1980’s where the volatility was driven by a fractional Brownian motion 
with a Hurst parameter, H £ (0.5,1) [13, 45]. Also, a variety of research is conducted recently on the 
fractional Black-Scholes models [12, 8]. In a recent contribution, the family of affine Volterra processes was 
introduced as solutions of certain stochastic convolution equations with affine coefficients which generalize 
the rough volatility models (see e.g. [38]). There exists also some efforts to enhance the flexibility of these 
models with a local volatility component [40]. 

The “roughness” of the RFSV models stems from the fact that the Hurst index of the fBm driver is less 
than 0.5 and so the sample paths of the log-volatility process behave more chaotically than the Brownian 
motion. These models have attracted a considerable attention recently and many researchers are now in¬ 
volved in modeling, pricing, hedging, portfolio optimization and Monte-Carlo simulation under this family 
of models (see e.g. [26, 18, 6, 27, 39, 56, 19, 36, 54, 42, 23, 24, 30]) as well as studying their asymptotic 
properties [21, 31, 33]. 

In a recent contribution, El Euch and Rosenbaum [20] found the characteristic function of the log-price 
increments in the rough Heston model using the connection between nearly unstable heavy-tailed Hawkes 
processes and fractional volatility models 1 . Si mi lar to the classical Heston model [35] in which the char¬ 
acteristic function is expressible as the solution of a Riccati differential equation, El Euch and Rosenbaum 
[20] showed that the characteristic function of the rough Heston model exhibits the same structure and could 
be expressed in terms of the solution of a fractional Riccati ordinary differential equation (FRODE). 

In order to apply the rough Heston model in practice, it is necessary to calibrate it to liquidly traded 
plain-vanilla European options by minimizing an error criteria defined as the square root of the differences 
between the model prices 2 and the corresponding market quotes (see e.g. [19] for more details). Since by 
numerical approximation of the solution of the above mentioned FRODE we have the characteristic function 
of the log-price at hand, the choice of a suitable nonlinear fractional ODE solver has a profound effect on 
the viability of the overall procedure in terms of speed, stability and accuracy. 

Among the existing numerical approaches for solving fractional differential equations, those based on 
finite differences [51, 15, 48, 16, 68], finite elements [70, 53] and spectral methods [64, 66, 61, 41] have 
gained much popularity in recent years. There exist other schemes to solve nonlinear fractional ODEs such 
as B-Spline collocation [41] which result in a system of nonlinear equations (see e.g. [50, 4] for an overview 
of numerical methods for general fractional ODEs). Due to the non-local and weakly singular nature of the 
fractional differential operator, the computational costs and storage requirements of numerical methods for 
FODEs are higher than their integer-order counterparts [50]. 

The fractional Adams method (see e.g. [15, 48, 16]) as a widely-used difference-based scheme is a 
predictor corrector type method based on Volterra integral equation formulation of the fractional ODE 

'Taking a different approach through the theory of finite-dimensional deterministic convolution equations, Abi Jaber. et al. [38] 
proved that, under specific conditions on the kernel function and the coefficients of the equation, the affine Volterra processes and 
their corresponding Riccati-Volterra equations have a unique global solution. In particular, they show that the rough Heston model 
has a unique solution and its characteristic function is determined in terms of the unique solution of its associated fractional Riccati 
equation. 

-Implied volatilities could also be used in the calibration process. 
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which is derived from applying a fractional Euler method in the predictor step followed by a fractional 
trapezoidal formula in its correction step. El Euch and Rosenbaum [20] and El Euch, et al. [19] employed 
this scheme in their numerical experiments concerning the computation of the at-the-money skew under 
the rough Heston model. This method while being easily implementable and having a super-linear rate of 
convergence, appeal's not to be a good candidate to solve FRODEs appearing in the rough Heston model 
(see Section 3). Recent research by Gatheral and Radoicic [30] shows also that in the Adams method, “... 
a large number of steps and computation time is required to achieve satisfactory accuracy”. They conclude 
that “... such numerical solutions are inevitably much slower to compute than the closed-form classical 
Heston expression.” 

The failure of the Adams method is due mainly to the fact that the solution of this FRODE depends on 
a critical parameter and the Adams method fails to produce reliable results when applied to the problem 
with large parameter values due to a sharp transient behavior of the solution in this case. The problem could 
be partially resolved by adjusting the step-size of the base grid to the parameter used. However, this will 
result in a very fine grid for some parameter values and this in turn will slow-down the process considerably. 
Moreover, during the calibration procedure, we don’t have any control on the intermediate parameter values 
that the optimization algorithm produces and there is no straightforward formula which relates the proper 
step-size to the model parameters. 

In order to circumvent the above mentioned limitations of the fractional Adams method, Gatheral and 
Radoicic [30] proposed a rational approximation of the fractional Riccati solution, based on Pade approxi- 
mants and valid in a region of its domain relevant to option valuation. Callegaro et al. [32] also studied the 
radius of convergence of a specific short-time series expansion and applied Richardson-Romberg extrapo¬ 
lation (outside the convergence domain) to obtain a “promising and encouraging” hybrid numerical scheme 
for all times. 

In this paper, we propose an alternative computational scheme for solving FODEs based on iterative 
quasi-linearization (also called Newton-Kantorovich iteration) followed by spectral collocation to solve the 
resulting linear FODE at each iteration using fractional (non-polynomial) Jacobipoly-fractonomials. Quasi¬ 
linearization is a powerful nonperturbative approximation technique to solve highly nonlinear ordinary and 
partial differential equations which produces a rapidly convergent sequence of iterations (see e.g [7,44]). On 
the other hand, (orthogonal) poly-fractonomials as solutions to (regular) fractional Sturm-Liouville eigen- 
problems are natural candidates to be used as basis functions for expanding the solution of FODEs. It is a 
well-known fact that spectral methods are “infinite-order accurate” if the expansion functions are properly 
chosen [22]. The complete theory of these bases could be found in [62] and other related works (see e.g. 
[63, 64, 67, 69]). 

The proposed scheme combines the generality of quasi-linearization with the accuracy of spectral col¬ 
location to approximate the fractional derivative operator of any order using a global spectral expansion. 
Other advantages of the proposed numerical scheme could be summarized as follows: 

(i) avoiding the numerical solution of systems of nonlinear equations; 

(ii) allowing the calculation of collocation matrices for the solution and its fractional derivative only once 
inside the calibration process and the Newton iterations; 

(iii) producing well-conditioned matrices; 

(iv) availability of the numerical analysis for the scheme. 

In order to ensure the validity of the proposed method and to guarantee the specified convergence rate, 
we first prove some regularity results on the linearized problem arising from the FRODE and then employ 
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the proposed scheme to solve a practical calibration problem from the SPX options market in which we 
use the implied volatility of the index options to measure the accuracy of fitting process. The obtained 
results confirm the efficiency and stability of the pricing engine which is based on the new FRODE solver. 
Comparisons with the Adams method show that the time and memory requirements scale suitably with 
increasing the number of option contracts. 

We organize the rest of the paper as follow. In Section 2, the rough Heston model is briefly reviewed. 
Section 3 discusses the fractional Adams method and its shortcomings in solving the fractional Riccati 
differential equation appealing in the characteristic function of the log-price in the rough Heston model. 
The alternative numerical method to solve this problem and its error analysis arc introduced and discussed 
in Sections 4 and 5. Finally, as an illustration of the efficiency of the method, option pricing and calibration 
of the rough Heston model is presented in Section 5. We conclude the paper with some final remarks and 
open questions for future research in Section 7. 


2. The Rough Heston Model 


We arc concerned in the remainder with the rough fractional Heston model describing the coupled 
dynamics of the underlying asset, S, and the volatility, V, by 

dS, = S,VV,dWt, S 0 = 5o€M+, 

V, = Vo + J-^ f , (t-s) a - l X(Q-V s )ds+=j-^ f {t-s) a ~ i 'k\yJV s dB s , V 0 eM + , 

1 (a) jo l (a) Jo 

in which A,, 0, Vo and v are some positive constants, { W ,; J > 0} and {B,:t > 0} are two correlated Brownian 
motions with (dW t ,dB t ) = p dt and a G (1 /2,1) is a parameter determining the smoothness of the volatility 
sample paths (see [20] for more details). This model benefits from all the advantages of the classical Heston 
model while being interpretable in terms of the micro-structural foundations of the market 3 . 

It is shown (see e.g. [20]) that, for — l/v/2 < p < 0, the characteristic function of the log-price, X t = 
log(S,/S 0 ), is given by 


\| f(a,t) •= K[e ,<!X '] = exp 


Q)il l h(a,t) + Vq/ 1 a h(a,t) 


Mt >0, 


( 1 ) 


in which h(a,t ) solves the following fractional Riccati differential equation: 

l n v ) 2 

D a h(a,t) = -(— a 2 — ia) +X(/npv — \)h(a,t) + ^ ^ h 2 (a,t), f€[0,r], (2) 

7 1_a A(a,0) = 0. (3) 


In (2) and (3), D r and l r denote, respectively, the Riemann-Liouville fractional order derivative and integral 
operators (see e.g. [4]) of order r € (0,1), defined by 

£>7(0 = r(hr)J,l (, ~ srrfis)ds ’ 

= J-T [ (t-Sy~ 1 f(s)ds. 
r(r) Jo 


3 More precisely, assuming a high degree of endogeneity of market structure, absence of arbitrage opportunities, buying/selling 
asymmetry and presence of meta-orders as the main futures of the market micro-structure, El Euch et al. [18] arrived at the rough 
Heston model which exhibits both leverage effect and rough volatility. 
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The existence and uniqueness of a continuous solution to (2)-(3) is studied in [20] and [38]. 

The FRODE (2)-(3) has two main features. First, it is a nonlinear complex-valued differential equation 
involving a fractional-order derivative operator with no known analytic closed-form solution. Second and 
most important, the parameter a in the equation has a profound effect on the solution behavior (see Figure 
1). In fact, in order to calculate the option prices using the characteristic function, we need to find the 
solution numerically for values of a with large absolute values (see Section 6 for more details). 


Remark 1. In order to solve the fractional Riccati differential equation (2)-(3), it is possible to work with 
its equivalent integral equation formulation of the form 


h(a,t ) 


l 1 ~ a h(a,0) 

T(a )H~ a 


1 

+ m) 



F(a,h(a 1 s))ds , 


(4) 


in which 

1 , ( h /) 2 , 

F(a,u ) = - (—a — ia) + X(iapv — 1 )u+ u. (5) 

Note that based on condition (3), we could easily conclude that h(a, 0) = 0. Equation (4) is a parametric 
second-kind nonlinear weakly singular Volterra integral equation. 


Remark 2. Since the initial condition of the FRODE (2)-(3) is not given explicitly in terms of the function 
itself (but in terms of its Riemann-Liouville fractional integral), we could write this equation with the same 
right hand side but now in the language of Caputo’s fractional derivative [14], defined by 


D *«(0 = 


1 


r(i-a) 


f ( t-s)- 
J 0 


t du(s) 
ds 


ds , 


and satisfying the relation 


Dy(t) = D a u(t) - r J_ a ^ (°). 

Since the initial condition of the fractional differential equation in the Caputo sense is given explicitly in 
terms of the function itself, it is usually easier to develop numerical schemes for fractional differential 
equations in the Caputo sense. However, it is easily seen that when the initial conditions arc zero, the two 
problems arc equivalent and so numerical as well as theoretical analysis provided for the Caputo sense frac¬ 
tional equation are also valid for the equation in the Riemann-Liouville sense under zero initial condition. 
Based on Remark 1 which implies the zero initial condition, from now on and for ease of notation, we use 
D a interchangeably for both Caputo and Riemann-Liouville cases. 

Remark 3. Note that if h(a,t ) solves the integral equation (4) for the function F defined in (5) and 
h(a, 0) = 0, then h{a,t) := h(—ia,t ) solves the new integral equation 

h(a,t) = / (t — s) a ~ 1 F(a,h(a,s))ds, (6) 

r(a) Jo 


with 


F(a,u ) 



a) + X(apv — 1 )u + 



(V) 


Based on Remark 3 which guarantees the transferability of results from one case to another, we present 
the numerical results in Section 3 based on the complex-valued problem (2)-(3) and the theoretical error 
analysis in Sections 4 and 5 based on the real-valued version (6)-(7). 
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3. Fractional Adams Method and its Limitations 


In this section, we first briefly review the essential ideas of the fractional Adams method and then 
illustrate some of its limitations in solving the fractional differential equations mentioned in the previous 
section. This method is one of the efficient schemes to solve equation (2)-(3) and is a generalization of the 
one-step Adams-Bashforth-Moulton algorithm to fractional-order integral equations. Since the function F 
in (5) is nonlinear in its second argument, we apply a predictor-corrector algorithm to solve the problem. 
Detailed derivation of this method is explained in [16] and [15]. 

In order to briefly outline this method for solving the equation (4), consider a subdivision of the time 
interval [0, T] into N sub-intervals of the form [ 4 , 4 + 1 ] i n which 4 = kA, k = 0, • • • ,N and A = T/N. To 
simplify the notation, we define g(t) = F(a,h(a,t)). Taking the initial condition (3) into account, we arc 
looking for a numerical approximation to the solution of the following equation: 

[ (t-s) a ~ 1 g(s)ds. (8) 

l (a) Jo 

The main idea is to use product trapezoidal quadrature formula to approximate the integral on the right- 
hand-side of ( 8 ). In other words, at any time step 4 + 1 , we apply the approximation 

rtk+i ftk +1 , 

/ (4+1-*)“ g(s)ds& / ( 4 + 1 - 5 )“ 'g k+ \{s)ds. (9) 

Jo Jo 

in which g k +\ is the piecewise linear interpolant of g with nodes and knots chosen at 4 , 0<j<k+l: 

8k+i{t) = tj+1 _[ g{tj) + 1 j 7 g{t j+ 1 ), t G [ 4 , 4 + 1 ), 0 < j < K. 

4+i 0 O+i 0 

After some lengthy but straightforward calculation, we arrive at 

rtk+i ~ n 

J o {t k+ i-s) a ~ l g k+l (s)ds = ^ [ Y, «M+1S(0) +g(0+i)] 

in which 

«o,jt+i = k a+l - (k-a){k + l) a , 

a j,k +1 = (k — j + 2)“ +1 + (k — j) a+l — 2(k — j +1)“ +1 , 1 < j < k. 

So, the collector formula can be written as 

^ A“ ^ ^ 

Had (a, 4+1) = r , F(a,h p (a 1 t k +i))+ Y a jJ+^ F ( a ^ h Ad(a, t j)) , 

1 [a + zj j_ 0 

in which the functional equation of the Gamma function is used twice (which yields r(a)a(a + 1) = 
r(a +2)). To finalize the scheme, we need the predictor formula to calculate h p (a,t k +i)- To do this, 
we approximate the integral term in left-hand-side of (9) using the product rectangle rule as 

r 1 1+ 1 _. A“ JL 

/ ( 4+1 - s) a g(s)ds w — £ b j}k+ ig(tj) 

jo u /=0 

in which 

bj,k+ 1 = (fc+1 -j) a ~(k-j) a , 0 <j<k. 
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This will result in a predictor approximation for h p (a,t k + 1 )) of the form 

^ A a k _ 

h p {a,t k+ 1 ) = r(a+1) Y, b j,k+lF(a,h Ad (a,tj)). 

Based on the error analysis presented in [16] (see also [48]), it could be shown that for values of a € 
[0.5,1), if we have h(a , •) € C 2 [0, T] for all a G R, then 

I h(a,fj)-h Ad (a,tj)\ < CtJ~ l A 2_a , 

and in particular: 

max 0 <j< n \h(a,tj) -h A d(a,tj)\ = 0(A), 


and 


maXt j€[e j]\h(a,tj)-h Ad {a,tj)\ = 0( A 2 “), 
for any £ € (0 ,T], 

However, theoretical as well as numerical evidence shows that the right-hand side derivative of the 
solution, h(a,t ) w.r.t. t as t — >• 0 + goes to infinity, when a —>■ °o (see Appendix A for a theoretical explanation 
as well as Fig. 1 for a numerical illustration). This explains why the fractional Adams method explodes and 
fails to produce valid results for large values of a and returns an undefined numerical result (see Table 1). 

The exploding instability of the fractional Adams method also stems from the fact that this scheme is 
conditionally stable (see Appendix B). As can be seen from Table 1, although refining the mesh can resolve 
the issue for some values of a , it just postpones the appearance of “NaN” to larger values of a. It is worth 
noting that this critical "large” value of a depends on the model parameters as illustrated in Fig. 2. In this 
figure, the cross section 9 l[h(a, 1)] is plotted for a = 0.6 and different values of the parameters p, X, and v. 
Wherever some parts of the plot are missing, the algorithm has returned undefined numerical results. We 
must note that the rightmost column of Table 1 shows the mean relative error (MRE) calculated by the 
relation 


MRE (A!) 


1 y-, \h(a. tj'j h A(i (a, tj ) | 
[ \h(a,tj) | 


( 10 ) 


for the fractional Adams method. Since the closed-form analytical solution of the FRODE does not exist, 
we have considered an approximate solution computed by the fractional Adams method with N = 12800 
time steps as a proxy for the “exact” solution. 

Based on the shortcomings of the fractional Adams method outlined above, it is natural to think about 
other alternative approaches to deal with FRODEs arising from rough fractional stochastic volatility models. 
As a first remedy, it is possible to use non-uniform grid points in the finite difference discretization using 
optimal mesh-sizes obtained from stability conditions of the underlying scheme. This topic deserves a 
separate attention in its own right and will be pursued elsewhere. The other alternative which is pursued 
here is based on the theory of quasi-linearization and spectral collocation which will be outlined in the 
remainder. 
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5?[/i(a, t)] 


a = 0.6 a = 0.8 a = 1 



Figure 1: The solution of fractional Riccati equation for p = —0.5, V = 0.5, X = 0.2 and different values of 
a. 

Table 1: Approximate solution of the fractional Riccati equation (2)-(3) at T = 2 for parameter values 
p = —0.5, V = 0.5, X = 0.2 a = 0.6 and different values of a using the fractional Adams method. N At j is 
the number of time steps. 


a 

NAd 

X{h Ad (a,T)} 

3{h A d(a,T)] 

time (s) 

MRE(A Af/ ) 

10 

50 

-49.0659150695 

14.8431562460 

0.03 

0.3806 


100 

-49.0744759134 

14.8415100884 

0.04 

0.1340 


200 

-49.0772023771 

14.8410382227 

0.12 

0.0486 


400 

-49.0780825041 

14.8408958427 

0.33 

0.0183 


800 

-49.0783690070 

14.8408514172 

1.17 

0.0070 


1600 

-49.0784627563 

14.8408372543 

4.40 

0.0026 


3200 

-49.0784935316 

14.8408326774 

17.09 

8.3e-04 


6400 

-49.0785036544 

14.8408311858 

72.20 

3.9e-07 


12800 

-49.0785069884 

14.8408306972 

349.75 

— 

100 

50 

-8.21733370e+02 

5.01209488e+02 

0.02 

0.3150 


100 

-8.16142423e+02 

4.63988366e+02 

0.04 

0.1042 


200 

-8.16515723e+02 

4.64180004e+02 

0.08 

0.0376 


400 

-8.16584730e+02 

4.64213767e+02 

0.36 

0.0138 


800 

-8.16601994e+02 

4.642224210e+02 

1.23 

0.0051 


1600 

-8.16606859e+02 

4.64224928e+02 

4.38 

0.0018 


3200 

-8.16608317e+02 

4.64225697e+02 

17.54 

5.6e-04 


6400 

-8.16608769e+02 

4.64225939e+02 

78.16 

3.7e-06 


12800 

-8.16608913e+02 

4.64226017e+02 

327.96 

— 

1000 

1600 

NaN 

NaN 

_ 

_ 


3200 

-8.61056342e+03 

4.96535070e+03 

18.27 

5.2e-04 


6400 

-8.61052875e+03 

4.96539020e+03 

59.32 

5.5e-05 


12800 

-8.6105321 le+03 

4.96539270e+03 

233.02 

— 

1500 

3200 

NaN 

NaN 

_ 

_ 


6400 

-1.29406708e+04 

7.46541212e+03 

58.22 

1.8e-04 


12800 

-1.2940663 le+04 

7.46543919e+03 

230.09 

— 
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Figure 2: The fractional Adams method fails to produce valid results for large values of a for different 
parameter sets. Wherever some parts of the plot are missing, the algorithm has returned undefined numerical 
results. 


4. Quasi-Linearization 


In order to deal with the nonlinearity present in (2)-(3), we apply the iterative quasi-linearization 
(Newton-Kantorovich) scheme resulting in a linear fractional ordinary differential equation (FODE) at each 
step and then approximate the solution of this linearized problem using the Jacobi poly-fractonomials to 
cope with with the fractional derivative operator. 

Let us consider the nonlinear fractional ordinary differential equation of the form 

D a u(t) = /(?,«(/)), 0 <t<T, (11) 

k(0) = u°, (12) 


in which the function / E C((0, T] x M)) satisfies the following two assumptions: 

Assumption 1. For some 0 < a < a < 1, the two functions t a f(t,.x) and t a ^/(tjx) belong to C([0, T] x M)). 


Assumption 2. For some 0 < a < a < 1 and some L > 0, we have 
for all t € [0, T] and X] ,x 2 E M. 


< L\x\ — x 2 \ 


Remark 4. It is straightforward to verify Assumptions 1 and 2 for the fractional quadratic Riccati differ¬ 
ential equation with f(t,x ) = cq + c ],r + C2A' 2 for some constants co,ci and C2 . Obviously, for any a > 0, 
t a f(t,x ) and t a ^~ x f(t,x) = t a (c\ +2c2.r) belong to C([0, T] x M) and so the Assumption 1 is satisfied. More¬ 
over, we have 


Yx f(t M )-gf( tl x 2 ) 


< T a 




< 2c 2 T a \xi —x 2 \ 


for all t E [0, T] and xi,x 2 € M and so the Assumption 2 is also satisfied with L = 2c 2 T a . 
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Based on the fact that the function / is nonlinear in its second argument, quasi-linearization is proposed 
here as a suitable method to deal with the nonlinearity issue. This method introduces an iterative procedure 
which, starting from an initial guess, improves the solution at each iteration using the first term of the Taylor 
expansion of the underlying operator (see [2] for more details). In other words, the quasilinear iterative 
scheme is obtained by linearizing the equation around the n- th approximation by expanding the function / 
in a power series and retaining only the linear terms. In the case of fractional differential equation (11)-(12), 
it is defined by the recurrence relation 

= /(f,M„-t) + /2(f, - M,-t), W„(0) = K°, H = 0,1,2, - - - , (13) 

in which w„_i(-) and u„(■) arc successive approximations to the solution and fi{t,x) '■= We 

emphasize that for every n > 0, u n i (•) is known at the n-th iteration, either from the initial guess (n = 1) 
or from the previous steps (n > 1). Consequently, f(t,u n - 1 ) and fi(t.u n -\) arc single-variable functions of 
time, t. Therefore, despite the presence of / and fi, equation (13) is linear w.r.t. u n (■). The procedure is 
summarized in Algorithm 1. 


Algorithm 1 Quasi-linearization Algorithm 
1: Pick an initial guess uq and an error tolerance e; 

2: Set n = 0 and v = 1; 

3 : while v > E do 
4: Increase n by 1; 

5: Solve the linearized problem (13): 

(D a u n -1 )I)Un = f{t , u n - 1) - f 2 (t , Un -1 )u n -\ \ 

6: Update v by v := \\u n — u n -i || ; 

7: end while 


In the next subsection, we provide sufficient conditions on the initial guess, uo(-), under which the 
Algorithm 1 is convergent and the rate of convergence is also derived. Moreover, an example of a suitable 
initial guess is provided. 


4.1. Convergence Analysis of the Quasi-Linearized Iterative Process 

The proposed iterative procedure (13) could be considered as an operator Newton scheme to solve the 
weakly-singular integral equation 

m(Q = m°+ 1 / (: t-s) a ~ l f(s,uis))ds , (14) 

1 (a) Jo 

which is equivalent to finding the (unique) zero of the operator ‘J : C([0, T],M) —>• C([0,7’],M), defined by 


QFu){t) = u{t)-u°-f^-T f [t-s) a l f(s,u(s))ds. 

i (a) Jo 


(15) 


In order to guarantee the convergence of Newton iterations to the solution of (14), we need to show that 
the operator f is Frechet differentiable and it’s Frechet derivative, ‘f'f h) at the point u G C[0,T], defined 
implicitly by the relation 


lim 

\\ h\\-*0 


!f(u + h)-!F(u)-!Ft(h) 
M 


= 0 , 
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has a bounded inverse. It could be shown (see e.g. [59]) that the operator ‘J defined by (15) is Frechet 
differentiable everywhere and we have 

fiKt) = h (t) ~ fj-y [ (■t-s) a ~ 1 f 2 (s,u(s))h(s)ds , V/7 G C[0,F]. 

l (a) Jo 

In order to ensure the quadratic rate of convergence for the iteration (13), we need the following theorem 
(see e.g. Vijesh [59]) which provides sufficient conditions for the global convergence of the scheme. Before 
presenting the result, we review some features of the problem which will be used in the sequel. Choosing a 
suitable value for a satisfying Assumptions 1 and 2, we define the constants M,r\ G M by 

M= sup \t c f 2 (t,u 0 )\, (16) 

fe[o,r] 

and 


r| = sup \JFu 0 (t)\, (17) 

f€[0T] 

in which no G C[0, T] is the starting point of the iteration and the operator ‘J is given by 15. Consider now 
the following two assumptions: 

Assumption 3. The constant M defined in (16) satisfies the inequality 

T(1 -o)T a ~ a M 
r(l+oc-a) < ’ 

for some 0 < a < a < 1. 

Assumption 4. The constants M, r| and L defined respectively in (16), (17) and Remark 4 satisfy the 
inequality 

Lr( 1 + a - a)r(l - a)r“-°r| 1 

(r(l + a-a)-r(l-a)r a “°M) 2 - 4’ 

for some 0 < a < a < 1. 

Theorem 1 ([59], Theorem 3.1). Let u o(-) G C[0, T] and suppose that Assumptions (3) and (4) arc satisfied. 
Then the fractional differential equation (11)-( 12) has a unique solution u* in U (no, r), the closed ball with 
center no and radius r, with 

2T(1 + a-a)r| 

r ~ T(1 +a —a) -MY(\-a)T a ~ a ' 

Furthermore, the quasilinear iterative scheme (13) starting from no converges quadratically to u* and for 
each n G N, we have 

l|w*-n„|L<c 2 -"e 2 ' ! , 

in which 

Q _ 4Lr(l + a-a)r(l-a)r a -°r| ^ 

“ [T( 1 + a - a) - MT( 1 - a) r a ~°] 2 ~ 

In the above relations, C is some constant and WfW^ = sup f€ r 0 7 \f(t)\. 
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Proof. Vijesh [59] stated the above theorem for the Caputo fractional derivative operator. However, his 
proof only depends on the integral representation of the problem, which is the same for both Caputo and 
Riemann-Liouville cases, up to a fractional power function of the independent variable which vanishes un¬ 
der zero initial condition (see Remark 2 and [14] for more details on the integral representation of FODEs). 
So the result is also valid for the Riemann-Liouville fractional derivative operator with a zero initial condi¬ 
tion. 

□ 


Remark 5. (A Suitable Initial Starting Point) In the case of fractional quadratic Riccati differential 
equation of the form 

D a u = co + cin + C 2 n 2 =: f(t,u ), (18) 

with initial condition, u( 0) = 0, the quasilinear iterative method takes the form: 

D a u„ - (a + 2c2U„-i(t))u n =C 0 -C 2 [u n -i(t)] 2 , u„{ 0) = 0. (19) 

Let us denote by u( ■) the solution of corresponding (integer-order) ordinary Riccati differential equation 

u(t) = f(t,u{t)), 
u( 0) = 0, 


with the exact solution 

u(t ) = 


2c 0 


\Jc\-4c 0 c 2 coth (\Jc\-4c 0 c 2 t/2) -cj 
Now, if we define 


M = sup |/ 2 (/,w(f))|, 

(€[0,7’] 

rj = sup \!Fu(t)\, 

te[0,T] 


and select £1 and £2 satisfying 


0 < 
0 < 


r(i+a) 

MT a ’ 

(r(i+a)-Mr a ) 2 
8c 2 r(l+a)r a ri ’ 


then by setting £ = min{£i,£ 2 } and choosing lift) := £u(t), all of the conditions of Theorem 1 arc satisfied 
by a = 0, M = £M and r| = Erf. 


Based on the fact that the lineai’ized problem (13) does not admit a closed form exact solution in the case 
of fractional Riccati equation, we propose a spectral collocation method to approximate u n {-) numerically 
which will be discussed in the next section. 
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5. Spectral Collocation using Poly-Fractonomials 

Spectral collocation methods arc efficient and highly accurate numerical schemes to solve linear and 
nonlinear ordinary and partial differential equations (see e.g. [37] and the many references therein). They 
arc based on the expansion of the solution in terms of some globally supported basis functions in which 
the expansion coefficients arc obtained in such a way that the differential equation is satisfied exactly at a 
set of so-called collocation points. The fundamental unknowns arc the solution values at these points and 
the most appropriate choice of collocation points is the set of Gauss-Lobatto-Chebyshev quadrature nodes 
based on the theory presented in e.g. [34]. 

The popularity of these methods stems from several advantages which they have over common finite 
difference methods. First, they have the potential for rapidly convergent approximations. The dual repre¬ 
sentation in physical and transform space allows for automatic monitoring of the spectrum of the solution, 
providing thus a check on resolution. If certain symmetries exist in the solution, spectral methods allow 
the exploitation of these symmetries. Finally, the methods have low or non-existent phase and dissipation 
errors [37]. 

The linearized problem (13) could now be solved using spectral collocation methods. Note that this 
equation contains a global fractional derivative operator, suggesting the use of global numerical approx¬ 
imation schemes such as spectral methods, rather than the local ones such as fractional Adams method. 
But still, due to the presence of the fractional derivative operator, the choice of suitable and efficient spec¬ 
tral basis functions is crucial. Zayernouri and Karniadakis [63] have shown that Jacobi poly-fractonomials 
as analytical eigen-solutions of the regular fractional Sturm-Liouville problem (RFSLP) of first kind are 
suitable choices with favorable properties. In their consecutive work [64], they also have introduced the 
shifted Jacobi poly-fractonomials and their fractional derivatives. So, we are provided with a natural basis 
function 4 to solve our linear fractional Riccati differential equation using a spectral collocation scheme. 


5.1. Fractional Basis Functions 

The employed spectral collocation method relies on Jacobi poly-fractonomials explicitly defined as 


V»(x) = (x+1 yp n ^(x), xE [— 1 , 1 ], 

in which P n “'(‘{x) is the standard Jacobi polynomial with p G (0,1). Recall that for a,b > — 1, P k ’ b (-) denotes 
the Jacobi polynomial of degree k which has the explicit representation 

pU.b / \ _ y-' ( — 1)* " l (l+b)k(l+a + b)k+ m / 1+x\ 

k X J^ 0 m\(k-m)\(l+b) m (l+a + b) k \ 2 ) 


where for 9 € M, (0)o = 1 and (9),- = 9(9 + 1) • • • (0 + i— 1) for i = 1,2,. 

Applying the affine mapping x(t) = 2 t/T — 1 to transform [0 ,T] into [—1,1], gives us the shifted basis 
functions: 

n(t)=n( X {t)) = (lYpp^f(x(t)), t g [o,r], (20) 


Using the fractional Riemann-Liouville derivative of power functions (see e.g. [57]) of the form 

r(i+v) 


mt v 1 = 2 u 1 ' > f-u 
lJ r(i+v-p) ’ 


v > 1, 


4 The use of these basis functions will result in exponential error decay phenomena which is theoretically proven for general 
Fourier-based spectral or pseudo-spectral methods (see e.g. [22]). 
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the lineality of the the fractional operators and a little algebra, the fractional Riemann-Liouville derivative 
of the transformed function (20) could be obtained as 

D'WM] = (f PD 

in which A ( - ) = P^' (■) denotes the standard Legendre polynomial of degree k. This formula states that the 
p-order fractional derivative of such poly-fractonomial basis functions of degree (n — 1 + q) arc multiples 
of the standard Legendre polynomials of integer degree (see [64, 65] for more details). 


5.2. Problem Discretization by Collocation 

The main motivation for the use of spectral methods in numerical analysis is the attractive approximation 
properties of orthogonal polynomial or (nonpolynomial) poly-fractonomial expansions. It is shown in [63] 
that Jacobi poly-fractonomials, {tP“(x)}“ =0 , arc orthogonal with dense linear span in I?_ a _ a [— 1,1], As 
will be shown in the proof of Theorem 3, the solution, u, of the fractional Riccati equation (18) and also 
the solution, u„, of the linearized problem (19), both belong to Zr a a [(). T] for each n and so they could be 
expanded in a shifted Jacobi poly-fractonomial series on [0, T] as 


«w = 

}=i 

oo 

U n {t) = 

7=1 


In the above expressions, the expansion coefficients arc given analytically by 

rT 


Sj = A 


rj-< \ 2a 1 


c7_“' a Jo 


u{t)(Q{t)$’j‘(t)dt, 


- ( 2 


rji \ 2a 1 


— / u n (t)(o(t)$j-(t)dt, 


CJ^i Jo 


( 22 ) 

(23) 


where the weight function is given by ( 0 (t) = (T —t) a t “ and the orthogonality constant, C ; a ct , is defined 

by 


^-a,-a _ II D -a,ct ||2 

C j ~ ll F j llL-l c 


(l-x) a (l+x) b P. 


-a,a 


(x)^ dx. 


It is well known that the magnitude of the expansion coefficients (22) and (23) will depend only on the 
degree of smoothness of the solution. In particular, if u and u n arc infinitely differentiable, then the q 7 and 
will decay faster than any polynomial of 1/ j. 

Among the many different ways to approximate the solution of (13) by a finite-term expansion, the 
spectral collocation method is based on finding an approximation u n jf(-) to u„(■) belonging to 


^ = {v\v = C%, tigPaKM)}, 


in which P^([a,b]) is the space of all algebraic polynomials of degree up to N defined on the interval [a.b]. 
Choosing the shifted Jacobi poly-fractonomials as a base for and setting p = a. the function u h .n(-) and 


14 



its fractional derivative could be expanded as 


AM-1 __ 

UnA*) = (24) 

7=1 

AM-1 _ 

D a [unAt)] = L Z,njD a [V?(t)], (25) 

7=1 

in which D“[iP“(7)] is given by equation (21). 

As our collocation points, we use the shifted Gauss-Lobatto-Chebyshev quadrature nodes in the interval 
[0,T] given by 

tj = —(xj + l), j = 1,2,— ,N + 1, (26) 

in which the jt/s are the standard Gauss-Lobatto-Chebyshev nodes defined by 
/ n 2 / — 1 \ 

ATy = C° s ( - T ) . J = 

Inserting the expansion (25) into the equation (13) and collocating the result at the nodes (26), the problem 
of finding the approximation, u h ,n{-) will t> e transformed into finding the coefficients q„ 7 which satisfy the 
linear system of equations of the form 

[D-F^\p]l n =b n ^. (27) 

In the above equation, we have 
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and then we employ the fractional spectral collocation scheme (27) to solve these linearized problems at 
each iteration. In the remainder of this section, we present an error estimate for the fractional spectral 
collocation scheme (27) and the complete numerical procedure discussed above. Based on the fact that the 
error in spectral collocation phase decays exponentially fast as the number of expansion terms increases 
(see e.g. [65]), we focus only on the accuracy of Newton quasilinearization as well as projection on the 
space of fractonomial basis functions. 


5. 3. Error Analysis 

In order to conduct convergence rate analysis of the approximation procedure based on the spectral 
expansion in (24), we use the following result from [17, Theorem 3.6] which provides an error estimate for 
the above spectral approximation method. 

Before presenting the result, we need to introduce some notation. Let 7 = [0, T] and define the function 
spaces 7?”' g (7) and ‘if „(7) respectively by 

B^(I) := {«!«(*> € L 2 nkMk (I),k = 0,1, • • • ,m}, 

■= {u\D S+l u € L] +Ull {l),l = 0,1 

in which id- denotes the fc-order derivative of u. Let also T^ g (7) denote the space of all functions u defined 
on the interval 7 with the corresponding norm |//|| / p ;/) < °o. The inner product and norm arc defined as 

( u U u 2)l 2 s (I) = J U[(x)u 2 (x)w y ’ 5 (x)dx, 

II II _ / n 1/2 

in which hT 8 (x) = (T —x)7x 8 and y, 8 > — 1. 

If we choose y= —a, 8 = a, m = 1 and s = 0 in Theorem 3.6 from Duan et al. [17], we obtain the 
following simplified version of their result for spectral projection. 

Theorem 2. Suppose for some positive integer q, that «(•) € c(7) nfl* (7) n ®V a (7), where a e 
(0,1) and define 

N +1 __ 

«»(•) = £ %&(■)■ 
i =i 


Then there exists a constant C independent of N such that for 1 < q < N + 1 we have 
||w - waHIoo < CN l ~ q \\D a+cl u\\ L 2^, 
and for q > N + 1, we have 


xr 1/2-a 

I u - u N |U < \\D a+l+N u\ 


12 

^N+l,N +1 


Proof. To prove the theorem, it is enough to let y = —a, 8 = a, m = 1, s = 0 in [17, Theorem 3.6]. 


□ 


Theorem 3 ensures us that the solution of the fractional Riccati equation, «(•), as well as the «-th Newton 
iteration, «„(•), satisfy the conditions of Theorem 2 for q = 2. We will use these two theorems in Theorems 
4 and 5 to obtain the convergence rate of our proposed scheme. 
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Theorem 3. For any positive integer n, if n„_i(-) G C(7), then the solution u n {-) of the quasi-lincari/cd 
problem (19) belongs to C(7) n7?l a _ a (7) n < B 7 L a _ a (I). Moreover, the solution u(-) of the fractional Riccati 
equation (18) with zero initial condition belongs to the same space. In particular, u G I?_ a _ a (7) and D a+1 u G 

Proof. We need to prove that at each iteration, u n {-) and its derivatives satisfy the following conditions: 

(i) ui k) eL 2 _ a+k _ a+k (i), £ = 0,l, 

(if) D«+V, €7^(7), / = 0,1,2. 


Note that we can rewrite equation (19) as the following Volterra integral equation: 
l 

Unit) = yy o (f-5) _(1 ^ a) ([ci +2c 2 u n -i{s)]u n {s) + co-C 2 U 2 n _ l {s))ds, n W (0) = 0. 

It is shown in [55] that u„ is weakly singular of order one (see Appendix A), which implies that u n G 
C[0,7’] n C 1 (0. T], Moreover, for every £ > 0, u\n\t) is absolutely continuous on £ < t < T and so is 
bounded on this interval. Similarly, being continuous on a closed interval, u n is bounded on [0, T]. Let us 
define 

Mq := sup u n (t), M\\= sup u\,\t). 
o <t<T e<t<T 


Moreover, u\,\t) = 0{t a 1 ), i.e. there are £i, A'i > 0 such that |ni^(7)| < K\t a *, V0 < t < £i. So we have 


r t-2a r (l_ a )2 


< °°, 


f Q (u n (s)) 2 (T — s)~ a s~ a ds < Mq {T-s)~ a s~ a ds = M 2 r(2 _ 2a) 
and so iin ] = u n G L 2 a _ a (7). Similarly 

f {ul, l \s)) 2 {T— s) 1 ^ a s 1 ~ a ds = f {(ui l \s)) 2 {T — s) 1 ~ a s 1 ~ a ds+ f {(ul?\s)) 2 (T — j) 1-< V~ a ds 
JO Jo J Si 

< K 2 r s «-HT- S y-^-«ds + M 2 f T (T-s) l - a s l ~ a ds 
JO J El 


1 2 —a 


+ 


2- 3+2a y/nT 3 - 2a r(2 - a) 

_ r(f-a) 

£i(£ir) 1_ “ 2 7q(2 —a,-l + a,3 —a;|-) 


—2 -|- oc 


oo 


where 2 / r i (a,b,s;z) is the hyper-geometric function. Therefore u„ 1 G 7^_ a+1 _ a+1 (7) and so u n G 7?l a _ a (7). 

In fact, from [55], we know that u^\t) = Kt a 1 + 0(t 2a ' 1 ). This means uf '(t) = ()(t a 2 ), i.e. there 
arc positive constants e 2 . 7f 2 such that \uk\t)\ < Kid 1 2 , V0 <t< £ 2 . Differentiating once and twice both 
sides of equation (19) with respect to t yields the expressions 

D 1+ % = 0(t a ~ l ), 

D 2+a u n = 0(t a ~ 2 ). 

Again, there are positive constants £i +oc ,7<fi +ce ,£ 2+a ,7f 2+a such that |D 1+a n„(f)| < K\ +a t a ~ l , V0 < t < £i +ce 
and \D 1+a u n (t)\ < K 2+a t a ~ 2 , V0 < t < £ 2+a . 


17 



Moreover, since u n belongs to C\0.T], so does D a u n . Now, the first and second derivatives of u n belong 
to C(0, T], so do D l+a u n and D 2+a u n . Denote 

M a := sup D a u n (t), M l+a := sup D l+a u n (t ), M 2+a ■= sup D 2+a u n (t). 

0<t<T £\+a.<t<T 


We have: 


[ (D a u n (s)) 2 ds <M 2 T < oo, 
Jo 


and so D a u n G Lq () = L . Moreover, 

r-T 
/ 0 

,2a 


J (_D 1+a w„(s)) 2 (T — s)sds = J (D ]+a u„(s)) 2 (T — s)sds + J (D l+a u„(s)') 2 (T — s)sds 


<K 


2 ef“ a (r + 2s 1+a (r-e 1+cc )) 


l+a" 


2a(l +2a) 


+M- 


r 3 -re 


l+a 


l+a 


T 3 —£ 3 
1 fc l+a 


< oo, 


and so D l+a u n G L\ j . Also we have 
rT 


r i ^ /‘^2+a o 

J (p 2+a u n (s)) {T— s) 2 s 2 ds = jI ( D 2+a u n (s )) ( T — s) 2 s 2 ds 

+ f {p 2+a u n {s)Y (T — s) 2 s 2 ds 


'£2+a 

✓ v'l c 2a 
^ A 2+a t 2+a 


^2+a 

1 T 2a 


, 1 

-r - + 


+ ^2+a 


p 5 Tp4 T^p 

c 2+a J c 2+a _ J c 2+a 


a & 2 +a — 2ae 2+ a 

-o3 y5 


+ 30 


< oo, 


and so D 2+a u n G 7+ 2 . So u„ G 21 2 a a , and the desired result follows. 

Following the same reasoning as above for the integral representation of equation (18), it is straightfor¬ 
ward to see that the solution, u, of the fractional Riccati equation (18) with zero initial condition belongs to 
C(7) n7?l a _ a (7) n $£ a _ a (7). In particular, u G L 2 _ a __ a (7) and D a+2 u G L 2 2 . □ 


It is worth noting in the passing that based on Theorem 2, the smoother the solution of a fractional 
differential equation is, the higher the rate of convergence of the spectral collocation approximation which 
will be obtained. 


Theorem 4. Stalling from the initial starting point introduced in Remark 5, the convergence rate of the 
proposed quasi-linearization spectral collocation method to solve the fractional Riccati equation (18) is 
given by: 

\\u-u n , N \L<C l 2- n Q 2n +C 2 N- 1 \\D a+2 u n \\i 2 , (28) 

11 11 2 

in which C\ and C 2 are some constants independent of n and N and 

8c 2 r(l + a)r ct r| 

[r(l+ a) -MT a } 2 

for M and r| as defined in Remark 5. 
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Proof. In order to establish the required result, we use the triangle inequality to write 
|| U U Hj N | | 33 ^ || U U n | | „ T" || U n ll n ,N I I oo > 
which by using conclusions of Theorems 1, 2 and 3, could be re-written as 

\\u-Uu,n\L <C l 2~ n Q 2n +C 2 N - 1 ||D a+2 M JL 2 , 

which proves the relation (28). □ 

Although we know that the Lr_ a _ a norm of D a+2 u n is finite for all values of n, it is not easy to find 
an explicit bound in terms of n for this term. However, if we compute the global error in the L 2 a _ a norm 
(Theorem 5 in the remainder), we could find an explicit bound in terms of n and N. 

The following Lemma provides some properties of Jacobi poly-fractonomials and the norms of weighted 
L p spaces which will be used in the proof of Theorem 5. 


Lemma 1. The following properties hold: 


(i) 


(ii) 


tpa 


(pa 

j 


_ ( 2 \ 2a_1 r -a,a 

~ \ T) S'-1 ’ 


' 2 


< tn 1 / 2 


rpa 

j 


(iii) \\u\\ L 2 < G3 1 / 2 | 


for all u G L 2 a _ a {1), 
r T . 


in which co(t) = (T — t) a t “ and 03 := / 0 ' a>(t)dt = - 


Proof, (i) Using the dehnition of Jacobi poly-fractonomials and C ; , , we have 

2 

—n 


ma 

j 


f) J Q ®(t)(t ap j°L a (x(t)j) dt 


2 \ 2a / T 


T J \2 

2 \ 2a-1 

T 


)/i( 1_x ) a o +s a (stm) dx 


p-a,a 

7-1 


2a-1 


c 


7-1 • 


The expressions (ii) and (iii) arc direct consequences of the Holder inequality. 


□ 


Theorem 5. Starting from the initial starting point introduced in Remark 5, the convergence rate of the 
proposed quasi-linearization spectral collocation method to solve the fractional Riccati equation (18) is 
given by: 

||w —Wnjv||r 2 <C 3 2-' 7 e 2 ' , (A + l)+C 4 iV- 1 ||D a+ 2 M|L , (29) 

5 L -a,-a 11 11 l 2.2 

in which C 3 and C 4 are some constants independent of n and N and 0 is dehned as in Theorem 4. 
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Proof. First, we note that 


\U ^n,N 11 jJL 


oo _ N+l _ 
7=1 7=1 


°° „ AM-1 __ 

E ^7+Effi-W)^ 

t 

AM-1 


r2 

^-cc,-a 


7 


< 05 1 / 2 ||M — Mjvlloo + 


< ||U u N\\ L 2_ a a + ^ I ^,j %n.j\ 

7=1 

0 \ ( 2 a— 1 )/ 27 V +1 

f E IWA/I( C JT) 1/2 

17 j-i 

/o \ (2a—1)/2 At+1 

< tn^CAr'IlD^n^^+^J i\^j-^j\(CjT) 1/2 , (30) 

where the last two lines are obtained using Lemma 1 and Theorem 2. Now, in order to bound the summation 
expression in the above inequality, note that for any j > 1 , we have: 


~Snj\ = 


j, \ 2a-1 | ,.j 


CJ^J o 


u(t) — u n (t))(o(t)¥j L (t)dt 


K 2 


rj-i \ 2a 1 -j^ 


c 


-a.a 

7-1 


u — u„ 


rpa 

7 


S 2 


2a-1 


' m \ Z.U 1 -i 

< -i-ca-e 2 ^ 1 / 2 


c 


-a.a 

7-1 


ipK 

7 


Lh 


a,-a 

/y\ 2a-l i / o \ (2a—1)/2 

(l) F^c,2-e^>/^ ^ (c-T)'/ 2 , 


- Uj 


(31) 


where the last two lines are obtained from Theorem 1 and Lemma 1, respectively. Putting inequalities (30) 
and (31) together, we get: 


\u-u n ,N \\ L 2 <®' /2 CN- l \\D a+2 u\\ L2 +03 1/2 Ci2"''e 2 "(iV+l), 


' 2,2 


which completes the proof. 


□ 


It is worth nothing that although the second term of the error bound is increasing w.r.t N, the whole term 
tends to zero as both n and N go to infinity, since the rate of decay of 2“”0 2 " is much faster than the growth 
rate of N + 1. Although the constants C 3 and C 4 in (29) are hard to estimate, it gives still information on 
how the error is related to the two basic parameters n and N appealing in the two phases of the algorithm. 
In order to achieve accuracy as well as efficiency, we could make a balance between these two sources of 
error contribution which could be done in at least two different ways: 

(i) determine n and N beforehand such that these two error terms are roughly of the same magnitude; 

(ii) first fix N and then determine n such that the error is smaller than a given tolerance level. 

We have pursued the second approach in this paper as will be made clear in the numerical experiments 
reported in the sequel. 
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Remark 6. The statement of Theorems 4 and 5 are also valid for more general families of nonlinear 
fractional Riccati differential equations with time-dependent coefficients of the form 

D a u(t) = Co(t) + Cl(t)u(t)+C2(t)u 2 (t), 

in which co(-), ci(-) and C 2 (-) arc twice continuously differentiable functions on [0.7']. These generalized 
FRODEs appeal - in modeling many different phenomena in fields such as physics, biology, physiology and 
optimal control [60, 5, 58, 1], This indicates that the domain of application for the proposed method is far 
beyond the rough Heston model. 

5.4. Numerical Experiments 

In order to show the effectiveness of the above proposed methodology, we solve the Riccati equation 
(2)-(3) this time with the spectral collocation and Newton-Kantorovich linearization. We have implemented 
the algorithm on a laptop PC with a Intel(R) Core(TM) i5-2410M CPU at 2.30 GHz and 6 GB RAM using 
MATLAB R2017b programming environment. The results are shown in Table 2. Here, Nf sc is the number 
of collocation points and si mi lar to Table 1, the rightmost column is calculated using the expression (10) 
with /?Ad replaced with hf sc and h(a.tj) is the approximation of the solution with Af sc = 160 collocation 
points. The reported times in the second last column of this table does not include the pre-allocation of 
the constants and calculation of poly-fractonomial basis functions. Also, the number of Newton iterations 
needed for convergence will be determined automatically within the algorithm which is reported in the third 
column of Table 2. Moreover, as illustrated in Figure 3, it properly approximates the solution for all values 
of a in combination with any other parameter sets and with any fixed number of basis functions without any 
explosion. Comparing the “time” and “MRE” columns of the two tables, we observe that the new method is 
very fast and at least quadratically convergent for small values of a and super-linearly convergent for larger 
values as illustrated in Figure 4 in which the reference line shows the linear rate of convergence. 

It is worth mentioning that the numerical results are in accordance with Theorem 1 which analytically 
guarantees that the convergence of the algorithm is independent of the value of parameter a. In the se¬ 
quel, we present the results of applying the above mentioned methodology in solving a practical financial 
engineering problem from SPX options markets. 


6. Option Pricing and Calibration 

In this section, we present some numerical experiments concerning the computational efficiency and 
complexity of the proposed method in pricing and calibration under the rough Heston model. For this 
puipose and in the pricing engine, we employ the Fourier option pricing technique of Fewis [47] and Carr 
and Madan [9]. 


6.1. Fourier Option Pricing Method 

Having the characteristic function, l P(a. T) in (1), we could employ the closed-form approximation to 
the pricing integral suggested by Fewis [46] (see also [47, 28]) which applies the inverse Fourier transform 
to provide an option valuation formula as a contour integral in the complex plane 5 . Based on this formula, 
the call option price could be calculated using the expression: 


C(S,K, T) = Se~ qT 


-VSKe~^ T / 2 [ 
n Jo 


91 [e iakx ¥(a — i/2,T)\ 

a 2 + \ 


da, 


(32) 


5 It should be noted that a similar approach to Lewis [46] is described by Lipton [52] in the context of FX option pricing. The 
Lewis-Lipton approach generalizes the method of Carr and Madan [9] since their dampening factor corresponds to the path of the 
contour integral in the Fourier transform. 
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Table 2: Approximate solution of the fractional Riccati equation at T = 2 and different values of a, and a = 
0.6, p = —0.5, X = 0.2, v = 0.5 using Jacobi poly-fractonomials and Newton iteration. Nf sc is the number 
of Jacobi poly-fractonomials used for the interpolation. The Newton iteration stops when max(9t(||w n — 
u n - 1 ||)) < 10 8 and n is the number of Newton iterations to achieve this precision. 


a 

-^fsc 

n 

X[hbc(a,T)] 

3[h isc (a,T)] 

time (s) 

MRE(Af sc ) 

10 

5 

4 

-49.0807046952 

14.8375090685 

0.003 

0.0031 


10 

4 

-49.0787853972 

14.8407277479 

0.004 

7.6e-04 


20 

4 

-49.0785452921 

14.8408126238 

0.005 

1.8e-04 


40 

4 

-49.0785130448 

14.8408280659 

0.006 

4.0e-05 


80 

4 

-49.0785091334 

14.8408301719 

0.01 

7.3e-06 


160 

4 

-49.0785086836 

14.8408304252 

0.03 

— 

100 

5 

7 

-8.04384294e+02 

4.62076565e+02 

0.006 

0.0564 


10 

6 

-8.17306929e+02 

4.63830399e+02 

0.006 

0.0154 


20 

6 

-8.16607059e+02 

4.64224853e+02 

0.007 

0.0028 


40 

6 

-8.16609239e+02 

4.64226616e+02 

0.008 

4.1e-04 


80 

6 

-8.16608991e+02 

4.64226076e+02 

0.01 

7.0e-05 


160 

6 

-8.16608982e+02 

4.64226055e+02 

0.04 

— 

1000 

5 

8 

-7.64716391e+03 

4.41402125e+03 

0.002 

0.2485 


10 

8 

-9.06210022e+03 

5.21767134e+03 

0.005 

0.1158 


20 

7 

-8.75517161e+03 

5.03733388e+03 

0.007 

0.0459 


40 

7 

-8.63261762e+03 

4.97061137e+03 

0.009 

0.0130 


80 

7 

-8.61126258e+03 

4.96493275e+03 

0.01 

0.0024 


160 

7 

-8.61053341e+03 

4.96539753e+03 

0.02 

— 

1500 

5 

9 

-1.13712037e+04 

6.56299912e+03 

0.001 

0.2881 


10 

8 

-1.37467576e+04 

7.92408707e+03 

0.002 

0.1378 


20 

8 

-1.32561485e+04 

7.63630767e+03 

0.005 

0.0593 


40 

7 

-1.30133620e+04 

7.49537785e+03 

0.007 

0.0200 


80 

7 

-1.29465085e+04 

7.46498434e+03 

0.03 

0.0044 


160 

7 

-1.29406868e+04 

7.4653601 le+03 

0.08 

— 


in which k = \og(S/K) + (r — q)T, and q is the dividend rate. 

6.1.1. Pricing Experiment 

In order to compare the pricing performance of the proposed method with the two available competitors 
in the literature (the fast Hybrid scheme of Callegaro et al. [32] and the fractional Adams method of El 
Euch and Rosenbaum [20]), we have prepared Table 3 based on the data presented in [32] for the other two 
methods. The details of parameter choices for the other two schemes is presented in [32] with enough detail 
but we only stress here that the step-size of the fractional Adams method is chosen such that the obtained 
results coincide (up to a tolerance level of 10 2 ) with the ones obtained from the fast hybrid scheme in term 
of the implied volatilities. A range of maturities from 1 day to 2 years and a wide spectrum of moneyness 6 
values is used in the comparison. The choice of N in FSC depends naturally on T (as shown in (28)) but we 


6 


The moneyness of an option is equal to the ratio between the strike price and the current price of the underlying asset. 


22 












Figure 3: The fractional spectral collocation method successfully generates the result for all values of a. 



Figure 4: Rate of convergence for fractional spectral collocation method. 


have not tried at all to do this choice optimally as it has a slight influence on the reported times and prices. 
The times reported in Table 3 for the FSC method are inclusive of all the computations and pre-allocations. 

Unless for the case of 1 day maturity (T = 1/256) and the strikes greater than 100%, the accuracy of 
the FSC and Hybrid schemes are comparable but the time efficiency of the FSC is better than the other 
schemes. For 1 day maturity and larger than 100% strikes, it seems that the FSC scheme has resulted in 
smaller prices for the call option which needs more careful investigation and analysis to identify causes and 
sources of this particular behavior. 
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Table 4: number of available option contracts corresponding to different strike prices considered in the 
calibration experiment for each maturity (in years). 


T 

0.1178 

0.2329 

0.6164 

1.3644 

Mj 

57 

5 

7 

3 


6.2. Calibration Experiment 

Equipped with a pricing algorithm, we tty to calibrate the parameters of the model using a suitable error 
criterion in terms of the model parameters. The calibration algorithm finds the parameters of the model for 
a set of market data by minimizing the sum of squared differences between the market implied volatility 
(IV) and those generated from the model. 

The option prices on SPX (S&P 500 index) are extracted from Chicago Board Options Exchange’s 
database (see www.cboe.com), on Feb 02 2017 at 13:39 ET. Only the options with positive volume arc 
considered in the calibration. The risk-free rate of return is set to be 0.02 and for each maturity, the dividend 
yield is calculated from the Put-Call parity formula: 

P(S,t) = C(S,t)-Sexp(-q(T -t)) +Xexp(-r(T -t)). 


From the six parameters of the model for each maturity T, two of them (0 and a) arc fixed to the 
values a = 0.6 and 0 = 0.1 by trial and error. We must note that if we don't prefix the parameter a during 
the calibration process, the time efficiency of the calibration procedure may be degraded slightly but the 
algorithm still performs well, without any explosion. The other four parameters ([p. v, Vo. AJ =: 0) arc now 
calibrated by solving the following optimization problem: 


Mj 

min V I 

i=i 


v® — v* l 

v T.i v T.i\ 


2 


„.-T <p<0 , 

v, Vo, X>0. 


in which Mj is the number of options with different strikes considered in our calibration experiment (see 
Table 4), V® = [V®,, V® 2 , ■ ■ ■ , V r %] is the vector of implied volatilities corresponding to option prices 
obtained from rough Heston model for the parameter set 0 and Vj = \Vj 1 , Vj 2 , • • • , Vj ;Vf; ] is the vector of 
market implied volatilities for maturity T. For each Maturity, the AARE and MARE columns in Tables 5 
and 6 arc the average and the maximum value of the relative error, across all strikes, defined respectively by 


AARE(T,0) 

MARE(T,0) 


1 ^|V®-V r *, 

m t h V h 


max 

i—l-Mj 




® _ V* I 
T,i v T,i\ 

V* 

v T,i 


Although the errors of the two methods does not indicate a notable difference, applying the fractional 
Adams method will cause some technical difficulties in the calibration process: to obtain the call price, the 
integral term appearing in the right-hand side of (32) should be truncated into a finite interval. For each 
maturity, the upper limit of this finite interval should be adjusted beforehand, otherwise the Adams method 
will explode for large values of a. Many trial and error experiments are conducted to achieve suitable 
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Table 5: Calibrated parameters 0 = [p,v, Vo, A,] using the combined quasi-linearization and spectral collo¬ 
cation. Parameters a = 0.6 and 9 = 0.1 arc prefixed. The initial guess for the optimization procedure is 
Po = —0.5 ,Vq = 0.5,Pq = mean(Vy) and Xq = 0.4. 


T 

fVfsc 

P 

V 

Vo 

X 

AARE 

MARE 

CPU time (s) 

0.1178 

5 

-0.6406 

0.6687 

0.0063 

0.5972 

0.0183 

0.0503 

251.0477 


10 

-0.6362 

0.4609 

1.0e-06 

0.9670 

0.0186 

0.0512 

313.2552 


20 

-0.6361 

0.4631 

1.0e-06 

0.9644 

0.0186 

0.0513 

469.4882 


40 

-0.6361 

0.4637 

1.0e-06 

0.9637 

0.0186 

0.0513 

1760.1473 


64 

-0.6361 

0.4638 

1.0e-06 

0.9635 

0.0186 

0.0513 

4098.1155 

0.2329 

5 

-0.7264 

0.6811 

0.0069 

0.6050 

0.0051 

0.0095 

12.3341 


10 

-0.7222 

0.6762 

0.0067 

0.5975 

0.0052 

0.0098 

14.6119 


20 

-0.7206 

0.6742 

0.0066 

0.5944 

0.0053 

0.0099 

27.2793 


40 

-0.7201 

0.6736 

0.0066 

0.5935 

0.0053 

0.0100 

75.4406 


64 

-0.7200 

0.6734 

0.0066 

0.5933 

0.0053 

0.0100 

205.5123 

0.6164 

5 

-0.7855 

0.6685 

0.0049 

0.6519 

0.0013 

0.0022 

11.9977 


10 

-0.7825 

0.6698 

0.0046 

0.6517 

0.0016 

0.0026 

15.3662 


20 

-0.7816 

0.6740 

0.0046 

0.6470 

0.0017 

0.0028 

29.8595 


40 

-0.7813 

0.6736 

0.0046 

0.6481 

0.0017 

0.0028 

87.4910 


64 

-0.7812 

0.6733 

0.0045 

0.6488 

0.0017 

0.0028 

240.4867 

1.3644 

5 

-0.8061 

0.7468 

0.0045 

0.7031 

3.3e-07 

3.9e-07 

11.3848 


10 

-0.8042 

0.7438 

0.0037 

0.6997 

1.0e-06 

1.6e-06 

13.7644 


20 

-0.8036 

0.7396 

0.0032 

0.7052 

3.7e-07 

5.1e-07 

25.0779 


40 

-0.8035 

0.7416 

0.0032 

0.7031 

5.3e-05 

8.3e-05 

64.2100 


64 

-0.8035 

0.7406 

0.0032 

0.7033 

3.8e-07 

4.8e-07 

183.4643 


Table6: Calibrated parameters 0 = [p. v. Vo, X] using the fractional Adams method. Parameters a = 0.6, 0 = 
0.1 arc prefixed. The initial guess for the optimization procedure is po = —0.5, Vo = 0.5, Vo = mean(Vf) 
and Xq = 0.4 for all T except for T = 0.6164 for which we have used po = —0.8. 


T 

<2 max 

N a 

N t 

P 

V 

Vo 

X 

AARE 

MARE 

CPU time (s) 

0.1178 

100 

512 

256 

-0.6387 

0.4537 

1.0e-06 

0.9827 

0.0166 

0.0538 

882.1240 

0.2329 

80 

256 

80 

-0.7182 

0.6955 

0.0072 

0.6161 

0.0015 

0.0025 

34.4670 

0.6164 

60 

256 

80 

-0.7831 

0.6925 

0.0057 

0.6533 

0.0011 

0.0021 

30.3321 

1.3644 

30 

256 

80 

-0.7997 

0.7639 

0.0037 

0.7586 

1.3e-06 

3.0e-06 

55.4553 
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(a) T = 0.1178 


(b) T = 0.2329 


(c)T = 0.6164 


(d) T = 1.3644 


Figure 5: Market and calibrated implied volatilities for different maturities using the fractional Adams 
method. 



(a) T = 0.1178 (b)T= 0.2329 (c) T = 0.6164 (d) T = 1.3644 


Figure 6: Market and calibrated implied volatilities for different maturities using the combined quasilin¬ 
earization and spectral collocation method. 


values of a max , N a and N, for different maturities, see Table 6. In contrast, the combined quasilinearization - 
spectral collocation method works properly well for all values of a which allows working with the integral 
without any formal truncation, using the “quadgk” function in the MATLAB programming environment. 

It must be noted also that we don’t need to use more than 5 collocation points, as it is not rewarding 
in terms of the performance (error) of the calibration procedure. In addition, the Adams method is time- 
consuming, since a “for” loop on the time grid is inevitable. As the numerical experiment shows, the 
proposed fractional spectral method reduces the CPU time at least to one third of the Adams method. 
This reduction will be significant for large data sets, as we see for T = 0.1178 where 57 options data are 
processed. It is worth noting that the fractional Adams method uses a partition consisting of only 512 or 256 
values of a in the integration procedure, while “quadgk” uses a partition of much bigger size. Moreover, 
collocation method seems to be independent of the number of poly-fractonomial basis, i.e. increasing the 
degree of Jacobi polynomials doesn’t improve the error and only five basis functions are enough to gain the 
same or better precision as that of the fractional Adams method. 

Finally, in order to show that the calibration procedure is robust with respect to the number of option 
contracts in the sample, we repeat the calculations corresponding to the maturity T = 0.1178 now with the 
aid of a sub-sample including 9 contracts selected evenly from the main sample. The results are reported 
in Table 7. They clearly demonstrate that the AARE and MARE error indicators are normalized deviation 
measures between the market-implied and model-implied volatilities. 

7. Conclusion 

In this paper, we have presented an efficient scheme based on quasi-linearization followed by spectral 
collocation using Jacobi poly-fractonomials to solve fractional Riccati differential equations. We provide 
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Table 7: Calibrated parameters 0 = [p,v, Vo, A,] using the combined quasi-linearization and spectral collo¬ 
cation. Parameters a = 0.6 and 9 = 0.1 arc prefixed. The initial guess for the optimization procedure is 
Po = —0.5, Vo = 0.5,Vo = mean(V^) and Xo = 0.4. The options data include a subsample of 9 contracts 
from the original 57 contracts. 


T 

Af S c 

P 

V 

Vo 

X 

AARE 

MARE 

CPU time (s) 

0.1178 

5 

-0.6309 

0.7539 

0.0073 

0.5601 

0.0205 

0.0315 

47.3779 


10 

-0.6263 

0.4844 

1.0e-06 

0.9902 

0.0208 

0.0332 

56.9425 


20 

-0.6262 

0.4867 

1.0e-06 

0.9876 

0.0208 

0.0333 

99.0243 


40 

-0.6262 

0.4873 

1.0e-06 

0.9868 

0.0208 

0.0333 

251.8572 


64 

-0.6262 

0.4875 

1.0e-06 

0.9866 

0.0208 

0.0333 

793.4031 


sufficient conditions under which the method is convergent and also obtain its convergence rate. We show 
the applicability of the method in pricing and calibration under rough Heston model recently introduced into 
the field of finance. We propose the new scheme as a reliable alternative to the fractional Adams method 
previously proposed to solve the problem. Comparisons with a hybrid scheme based on fractional power 
series expansions also shows the effectiveness of the proposed methodology. Among the possible ways 
to improve the method, it is possible to study the effect of adaptive selection of collocation nodes on the 
efficiency and stability of the method. 
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Appendix A. On Weak Singularity of Solutions to the Volterra Integral Equation (8) 

Using the results of Miller and Feldstein [55] on the regularity of solutions to Volterra integral equations, 
we show that h(a,t ) is weakly singular of order one, in the sense that 

= ( a2 ~ a ) 9 n C j-^ t a ~ l + °(t 2a ~~ 1 )- (33) 

at 2(1 +a) 

Definition 1 ([55], Definition 1), Suppose that v is a nonnegative integer and / is a function defined on 
(0, T] (or [0, T]). Then / is said to be weakly singular of order v, if and only if 

• /GC(0,r] ifv = 0or/ec v - 1 [0,r]nc v (0,r] ifv>0; 

• for each e > 0, the function f' y) (t) is absolutely continuous on the interval £ <t <T; and 
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the function a v defined by 


a v (t,f)=f(T) + f\f <y+1) (s)\ds, 0 < f < T, 

is of class L 1 (0, T). 

Theorem 6 ([55], Section 5). Suppose that x(t) is the solution of 

x(t) = fit) + f it — s)^~ p ^g(s,x(s))ds, 0 <t <T, 

Jo 

where v > 0 is an integer and 0 < p < 1. If / and g arc sufficiently smooth, then jc(-) is weakly singular of 
order V + 1 and tc( v+1 ) (t) = 0(t ~ p ) as t —>• 0. Moreover, 

= /|V+ ” 10) + Kt-r + 0(t''-<•). V> 1, 

in which K = 2g(0,/(0))(v + 1 - p) ■ ■ ■ (1 -p)/(v + 2-p). 

Using this theorem and setting v = 0, we get 

d ^ = Kr p + o{t l ~ 2p ), 

and so the relation (33) will be obtained easily. 


Appendix B. A General Stability Result for the Fractional Adams Method 

Let h be a fixed step-size and tj,j = 0,1,* • • ,n and w n j. j = 0,1, • • • ,n be the nodes and weights of a 
given quadrature rule for discretizing integrals. Stability analysis for a numerical scheme of the form 

n 

yn=gn + hY, w njK(t n ,tj,yj), n = p + 1, • • • ,N h , (34) 

7=0 

for a Volterra integral equation 

x(t)=f(t)+ [ K(t,s,x(s))ds , 0 <t<T. 

Jo 

could be conducted in several different veins. The most natural approach (called here the “perturbation 
stability analysis”) is concerned with analyzing whether or not a small perturbation in the initial values 
(yo,yt,-'' ■)’/)- i) and/or the right-hand side terms, g(t n ).n = 0,1,--- ,A), will cause a large error in the 
numerical solution. More precisely we have the following. 

Definition 2. Let the perturbations Ay/, j = 0, • • • ,p and A g n , n = 0. I. ■ • • , A), he admissible in the sense 
that 


A := sup Ay / + sup Ag n <°°. 

j=0A,-,p n=0,l,— ,Ni, 

Let also the values y„ + Ay„, n = p + 1, • • • ,2V/, satisfy the perturbed problem 


y„ + Ay„ =g n +Ag n + Y,w„jK(t n ,tj,yj + Ayj). 

7=0 
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Then the solution, y n of (34) is conditionally stable with respect to perturbations Ay/, j = 0,••• ,p and 
A g tll n = 0,1, • • • ,Nk, if and only if 

sup |Ay„| < CA, 

n=p+ 1,— ,N h 

for some constant C > 0. 

Based on Theorem 3.3 in [49], the fractional Adams method is conditionally stable in the sense of 
Definition 2. 

Appendix C. A Brief Outline of the Fast Hybrid Scheme of Callegaro et al. [32] 

In [32], the authors propose a novel hybrid scheme based on fractional power series expansion accom¬ 
panied with time discretization for solving (2)-(3). The main idea is to use an expansion of the form 

h(t)=h Xpy (t) = Y j a k t ka (35) 

k> 0 

and then applying the D a operator to this series expansion and inserting it into the equation. Doing so they 
obtain a recursive relation on the coefficients of the form 

“ 0 = 0 - “ l = f(T7T)' °»t. = (M 2 +m) r( X^ ) 1) . *>i. 

in which a* k 2 denotes the discrete convolution defined by 

k 

af = Yj a < a k-u k> 0, 

1=0 

They show that the convergence radius, R/, of the fractional power series (35) is bounded as 
0 < R { 

ower <Rh < ^uppen 

and the series is absolutely convergent for every t £ [0 ,Ri,). In order to control the error induced by truncat¬ 
ing the series expansion (35) at any order no, they derive some computable errors bounds and then propose 
a mixed (hybrid) scheme by merging two methods, one based on (35) for approximating h and the other 
based on time discretization to approximate P{h){t) and I l ~ a (h)(t). In the time discretization phase, they 
use the Richardson-Romberg extrapolation based on a conjecture on the existence of an expansion of the 
time discretization error. 
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